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Abstract. We numerically calculate the configurational entropy S con f of a binary 
mixture of hard spheres, by using a perturbed Hamiltonian method trapping the 
system inside a given state, which requires less assumptions than the previous methods 
[R.J. Speedy, Mol. Phys. 95, 169 (1998)]. We find that S con f is a decreasing function 
of packing fraction tp and extrapolates to zero at the Kauzmann packing fraction 
(f K ~ 0.62, suggesting the possibility of an ideal glass-transition for hard spheres 
system. Finally, the Adam-Gibbs relation is found to hold. 



1. Introduction 

The idea that the glass transition is driven by a decreasing of the number of accessible 
states upon lowering temperature (or raising density) is quite old [HEJE]. In this picture, 
if crystallization is avoided, an ideal glass transition is expected to happen at the point 
where the configurational entropy S con f (the logarithm of the number of states) vanishes. 
When the liquid enters in the the supercooled region, the dynamics becomes slower and 
slower and the particles get trapped for an increasingly longer time inside the "cages" 
made by their neighbors: the dynamics of the system can be successfully described 
as a "fast" motion of the representative point in the 3N configuration space inside 
metastable states, and a "slow" motion corresponding to jumps among states. Entering 
more in the supercooled region the number of accessible metastable states decreases and 
the extrapolation to zero of S con f defines the ideal glass transition. In experiments (or 
numerical simulations) the region close to the ideal glass transition is unreachable, due 
to the "apparent" arrest of the system at the so called glass-transition temperature (or 
density) when relaxation times become longer than experimental time scale. The above 
scenario has been shown to be valid for many interacting systems, based on smooth 
pair-potential (as Lennard- Jones liquids), for which the Potential Energy Landscape 
(PEL) approach pQElElCZllElElliniinilEllI^ and the replica method [H HE] have 
allowed to give numerical estimations of S con f and of the ideal glass transition. 

The overall picture is still not well established for Hard Spheres (HS), for which the 
existence of a glass transition is still an open question [T6l HZl UHl USl I2Q] . A particularly 
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important role seems to be played by the dimesionality of the system. In particular 
in d=2 dimension there are numerical [20l |2I] and theoretical [221 [22] evidences of the 
absence of a thermodynamic glass transition, while the opposite seems to be true in 
d=3 [TH [2l] . Moreover, the step- wise form of the interparticle potential does not allow 
a PEL analysis and different approaches have to be taken in consideration in order to 
calculate the configurational entropy S con j. In the past, attempts to estimate S con j 
have been performed based on different evaluations of the entropy in each single state 
[TBI 1251 [26] . Recently the replica method has been extended to the HS case for one- 
component systems [23| 124"]. 

In this paper we follow an approach, based on the Frenkel-Ladd method [27] and 
recently introduced in the study of Lennard- Jones systems [28J and attractive colloids 
[291 30J, to numerically estimate S con f for binary hard spheres. As in previous studies, 
the calculation is reduced to that of vibrational entropy S v ib, using the fact that the 
total entropy S can be decomposed into the sum of a configurational contribution S con f 
and a vibrational one S^: 

S — S con f + S v ib . (1) 

This expression is consistent with the idea that, at high enough density, there are two 
well separated time scales: a fast one, related to the motion inside a local state (the 
rattling in the cage), and a slow one associated to the exploration of different states. 
The total entropy S is obtained by thermodynamic integration, starting from the ideal 
gas state. The quantity S vi b is calculated using a perturbed Hamiltonian, adding to 
the original Hamiltonian an harmonic potential around a given reference configuration. 
Calculating the mean square displacement from the reference configuration and making 
an integration over the strength of the perturbation, it is possible to estimate 
the vibrational entropy [29]. The difference 5* — Smb provides an estimate of the 
configurational entropy S con f as a function of packing fraction (p (or density p). 

The main findings of the present work are the following, i) S con f is a decreasing 
function of packing fraction ip, and a suitable extrapolation to zero provides and estimate 
of the ideal phase transition point (Kauzmann packing fraction) ip K ~ 0.62. ii) The 
diffusivity D and configurational entropy Sconf are related through the Adam-Gibbs 
relation, in agreement with previous claims [IB] . 

2. The model 

The studied system is a binary mixture 50 — 50 of hard spheres, A and B, with diameters 
ratio o B j a A = \.2. The collision diameters are a AA =cr A , o BB =a B and o AB —{o ' A + a B )/2. 
The particles (iV=256) are enclosed in a cubic box with periodic boundary conditions. 
We use the following units: o B for length and rriA = tub — 1 for mass. Moreover we 
chose fee = 1 and h = 1. The density is measured by the packing fraction ip that is 
related to the number density p = N/V by ip=pTr(a A + cr^)/12. We analyzed systems in 
the range ^=0.425 — 0.580. Hard spheres systems depend only trivially on temperature 
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that sets an overall scale for the dynamics, consequently we perform all our simulations 
at T = 1. We performed standard event-driven molecular dynamics [31] and we stored 
several equilibrated configurations at different density. 

3. Diffusivity 

The diffusion coefficients D of the two species have been extracted from the long time 
limit of the mean square displacements (MSD) (r 2 (t))=N~ 1 ([r(£) — r(0)] 2 ) (r is the 
3N- vector of the coordinates): 

lim ~ QD . (2) 

To improve the statistical significance of the data, an average over 10 independent runs 
have been performed. In Fig [1] the mean squared displacements for the slowest cases, 
i.e. <p > 0.56, are presented for both species. It is clear that on increasing the density 
the MSDs develop the typical two step relaxation pattern. The first part of the MSD is 
purely ballistic while, at later stage, it reaches the diffusive regime, described by Eq. [2j 
Between this two regimes a plateau starts to develop. This is the clear indication of a 
caging effect. Each particle starts to feel the crowding of its neighbors and it is trapped 
in a cage for longer and longer time on increasing the density. The height of the plateau 
is the typical "cage diameter squared", 5r 2 . For both species we find 5r 2 = 3 • 10~ 2 cx 2 
for a = A or B represented by a dashed line in Fig. [TJ This is a clear evidence that the 
two species have the same caging effect. We shall return on the value of 5r 2 later on. 

In Fig. [2] the diffusivities D of A and B particles are plotted as a function of the 
packing fraction (p. Dashed lines in the figure are power-law fits D = C (ip c — y?) 7 
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Figure 2. Diffusivity for A and B particles as a function of packing fraction ip. The 
lines are power-law fits for ip > 0.53, C (ip c — tp) 1 , with: <,o c =0.583, 7=2.27, C=9.50 
for A particles; <y9 c =0.583, 7=2.47, C=11.66 for B particles. 



of the high packing-fraction data (ip > 0.53), as predicted by mode-coupling theory. 
The fitted parameters are: <y2 c =0.583, 7=2.27, C=9.50 for A particles and <^ c =0.583, 
7=2.47, C = 11.66 for B particles. We note that both diffusivities give rise to the same 
mode-coupling packing fraction <p c , in agreement with the prediction of the theory [32] 
and with previous simulations of the same model [33J . 



4. Configurational entropy 

We now turn to the calculation of configurational entropy. The method we follow to 
estimate S con f requires the computation of the total entropy S and vibrational entropy 
S V if,. The total entropy S is calculated via a thermodynamic integration from ideal gas 
and can be expressed as 

S(p) = S id (p)+S ex (p) , (3) 

where Sid is the entropy of the ideal gas and S ex is the excess entropy with respect to 
the ideal gas. For a binary mixture, the ideal gas entropy is: 

^ = ^-lnp-31nA + ln2, (4) 

where A = (2n(3h 2 /m)a is the De Broglie wavelength (h is the Planck's constant and has 
been set to unitary value), and the term In 2 takes into account the mixing contribution. 
The term S ex can be expressed in the following form 

Sex(p) = / % Pex , (5) 

1 Jo p z 
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Figure 3. Excess entropy S ex for the mixture of hard spheres as obtained from 
our simulations (symbols) compared with the analytical Carnahan-Starling expression 
(line) [351 E6]. 



with P ex the excess pressure. We extracted P ex from the zero density limit up to 
the densities of interest, performing numerical simulations and fitting the results of 
the pressure with a high order polynomial in p. In Fig. [3] we show the numerically 
calculated excess entropy S ex (symbols) together with the analytic estimate provided 
by the Carnahan-Starling (CS) equation of state, extended to hard sphere mixtures 
[3U [351 ESJ- We note that at high densities the CS equation of state overestimates the 
entropy of about 7%. This discrepancy, however, is not sufficiently significant to affect 
the resulting S con f values, in particularly closed to the glass transition. 

The method we use for the calculation of S v a is based on the investigation of a 
perturbed system 

(3H' = f3H + aN(r - r ) 2 , (6) 

where H is the unperturbed hard spheres Hamiltonian, a is the strength of the 
perturbation, r specifies the particles coordinates of a reference configuration and 
(r — ro) 2 = N" 1 YliLii^i — n>,i) 2 - The reference configuration ro is chosen from 
equilibrium configurations at the considered density (randomly extracted from the stored 
configurations obtained during molecular dynamics simulations). With this choice one 
is sure that the estimated vibrational entropy (see formula below) pertains to the correct 
state at the studied density. The vibrational entropy can be obtained from the formula 
(see Ref. (29] for details): 

W = L ^(r-r )> Q ,--ln(--) + -, (7) 

where a 0i oo are the lower/upper limit of integration, and (...) a i is the canonical average 
for a given a'. The choice of «o deserves few comments. If the system were confined 
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to move inside a given local free-energy minimum, for a correct estimation of S v n, one 
would take the lower limit a =0 in the integral in Eq. [71 As the system, at enough low 
value of a, begins to sample different states (the harmonic force due to the perturbation 
is no more able to constrain the system inside one state), cto has to be chosen in such a 
way that the system has not yet left the state: the underlying idea is that Eq. Ogives a 
correct estimation of S V ib until the system remains trapped in the state. An appropriate 
choice in our case seems to be a^=2 2 - 5 for all the densities , as, close to this point, one 
observes a crossover for all the investigated densities (more pronounced for low density 
data). In Fig. H] the quantity ((r — i"o) 2 ) a is reported as a function of a. An arrow 
indicates the chosen value ao=a^\ below which one observes the crossover associated 
to the exploration of different states. 

It is worth noting that different choices of «o are in principle possible, giving rise to 
different estimations of the vibrational entropy term. However, even though a kind of 
arbitrariness is present in the method, one can argue that a reasonable choice should be 
for values above (close to the crossover corresponding to the exploration of different 

(2) 

states) and below an upper value a at which one is sure that the system is still confined 
in a single state. The latter value can be estimated requiring that the MSD ((r — r ) 2 ) 
is always close/below the cage diameter squared 5r 2 ~ 3 ■ 10~ 2 o~ 2 (with a = A or B) 
(this has been estimated from the plateau of the mean square displacement, see Fig. [p. 

The chosen value in our case is a^=2 A - 5 (indicated by an arrow in Fig. H]). We then 
repeated the same calculation of S v n, using EqJTJ with the lower bound in the integral 

(2) 

a o — a o ■ m this way we obtain a lower and upper bound for the quantity of interest 
S V ib, by using respectively or in the expression of Smb in EqJTl 
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Figure 5. Configurational entropy S con f as function of packing fraction (p. Open 
symbols are data using ckq = 2 25 , full symbols using olq = 2 4 5 (see text). Dashed 
and dot-dashed lines are from Speedy [18] for binary and monatomic hard-spheres 
respectively. Thin full line is the analytical computation of Parisi and Zamponi for 
monatomic hard-spheres [24] . Thick lines are polynomial extrapolations of our data in 
the high packing fraction region, giving rise to the same Kauzmann packing fraction 
estimation for which S con f(ipK)—0- f K — 0.62. 



Fig. [5] shows S con f as a function of <p. The configurational entropy S con f is 
calculated using Eq. CD, where the two entropies 5* and S V ib are obtained from Eq.s 
[3] and [7J respectively. Due to the fact that the correct integral for the estimation of S v n, 
should be done from «o=0, but with the system always inside a given state, we have 
added to the expression in Eq. [7] the term a ((r — r ) 2 ) Q0 , corresponding to assume 
a constant value of ((r — ro) 2 ) below cto and using a zero value for the lower limit of 
the integral in Eq. [7J Fig. [5] shows the two estimates of S con f, corresponding to the 
two different values of a : a^=2 2 - 5 (open symbols) and a < Q^=2 4 - 5 (full symbols). One 
observes that the discrepancy between the two estimations decreases by increasing the 
packing fraction, suggesting that, at high density, the method used to calculate S con f 
is less affected by the choice of the free parameters entering in its evaluation. This is 
probably due to the fact that increasing the density the system tends to be more trapped 
in a local free-energy minimum. Indeed, it is only at high density that the method is 
expected to work better, due to the better definition of two time scales corresponding 
to local-fast and global-slow dynamics (see Fig. 1). At low density, instead, the two 
are less separated and this corresponds to a difficulty in the extrapolation for a — ► 
of the quantity reported in Fig. 4. The low density data show a more clear crossover 
on lowering a, and then a worst definition of state in this limit. As we are interested 
in the high packing fraction extrapolation, this fact do not affect our prediction on 
the Kauzmann density value. Also reported in the figure are the curves obtained by 
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Figure 6. The Adam-Gibbs relation D=D 00 exp[-NA/(TS con f)} (T=l) for the two 
species A and B: 1500=3.66, A=9.8 for A particles; Z?oo=4.59, A=10.5 for B particles. 

(2) 

The data of S con f are obtained with ao=aQ . 



Speedy [18J using a different method (assuming a particular form of the vibrational 
entropy, a Gaussian distribution of states and involving some free parameters) for the 
estimation of S con f, for monatomic (dot-dashed line) and binary (dashed line) hard- 
spheres (with the same diameters ratio 1.2 and composition 50 : 50 as in our case). It is 
worth noting that our method improves on Speedy's one, as, even though requiring some 
accuracy in the choice of the a parameter, has the advantage to be less affected by the 
presence of many free parameters and particular assumptions. We note that the data 

(2) 

of Speedy for the binary case do agree very well with our data with otQ=a Q , suggesting 
the possibility that the choice of ctQ=a i s m ore accurate for the estimation of S vi b 
and so of S con f. As a comparison, in Fig. [5] is also reported an analytic estimation of 
Sconf recently provided by Parisi and Zamponi [21] for monatomic hard-spheres. From 
the ^-dependence of the configurational entropy one can determine the packing fraction 
at which S con f extrapolates to zero, corresponding to the ideal phase transition point 
(Kauzmann packing fraction ip K ) S con f((p K )=0. Using a polynomial extrapolation [37] 
for the two set of data (corresponding to the different cto values) we obtain an estimated 
Kauzmann packing fraction value (p K ~ 0.62 (see Figj5]). It is worth noting that, even 
though the two curves are quite different, the estimated value of ip K is the same, again 
suggesting the robustness of the method in the high density region and then in the 
estimation of the Kauzmann packing fraction. 
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5. Adam-Gibbs relation 

In this Section we explore the validity of the Adam-Gibbs relation, linking dynamic 
quantities, like diffusivity, to S con f. In Fig. [6] we report the diffusivities D for A and B 
particles vs. the quantity {T Sconf)' 1 , with S con f obtained for the value ao=a^ . We 
find that the Adam-Gibbs relation 

D = D^e-^ 3 ^ (8) 

is verified (lines in the figure), with: 0^=3.66, A=9.8 for A particles; -Doo=4.59, A=10.8 
for B particles. A similar behavior is obtained using S con f calculated with a =a^ (not 
shown in the figure), with the values: -Doo=24.5, A=9.5 for A particles; -Doo=37.5, 
A=10.5 for B particles, suggesting that, in this range of diffusivity values, the AG 
expression is not able to discriminate between the two different estimations of S con f. 

6. Conclusions 

In conclusion, we have calculated S con f for binary mixture hard spheres, by numerically 
estimating the total entropy S (via thermodynamic integration from ideal gas) and the 
vibrational entropy S vi b using a numerical procedure based on Frenkel-Ladd method 
and recently applied in the analysis of Lennard- Jones systems and attractive colloids: 
the system is constrained inside a given "state" through an harmonic perturbed term in 
the Hamiltonian. We found, in agreement with analytical and simulation results in the 
literature, that S con f is a decreasing function of the packing fraction if, suggesting the 
possibility of a vanishing of S con f around the Kauzmann point (f K =0.62. Moreover, by 
studying the relationship between Sconf and the diffusion constant D, the Adam-Gibbs 
relation is found to reasonably hold for the analyzed system. 
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